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SUPERRESOLUTION  TECHNIQUES  AND  ISAR  IMAGING 


1 .  INTRODUCTION 

Radar  imaging  of  rotating  objects  has  been  an  area  of  interest  for  more  than  a 
decade,  and  several  descriptions  of  the  techniques/systems  involved  are 
available  in  the  literature  [1-5].  Essentially,  such  a  system  consists  of  a 
high-quality,  coherent,  pulse  dopoler  radar  feeding  its  received  signals  into 
a  suitable  range-Doppler  processor.  If  digital  processing  is  desired,  the 
radar  receiver  output  is  down-converted  in  a  synchronous  detector  to  baseband 
video  I+Q  (in-phase  and  quadrature)  signals  and  sampled  with  A/D  (analogue-to- 
digital)  converters.  Digital  range-Doppler  processing  may  then  be  carried  out 
in  convenient  discrete  Fourier  transform  (DFT )  devices.  The  resolution 
capabilities  of  these  systems  are  based  upon  the  well-known  principles  of  the 
conventional  windowed  Fourier  transform  both  in  the  Doppler  domain  and  the 
range  domain.  The  rotation  of  the  object  field  relative  to  the  radar 
generates  a  Doppler  frequency  gradient  and  permits  the  extraction  of  cross¬ 
range  resolution  that  is  much  finer  than  that  obtainable  by  the  radar’s 
antenna  beamwidth. 

This  paper  describes  an  investigation  wherein  high-resolution  optimal 
estimation  techniques  are  applied  to  the  radar  imaging  problem  area.  These 
methods,  sometimes  referred  to  as  "Superresolution"  techniques  [6],  are  a 
current  technology  that  is  of  considerable  interest  because  of  resolution 
capabilities  that  extend  well  beyond  the  conventional  windowed  Fourier 
transform.  They  have  been  sucessfully  applied  in  spectral  analysis,  seismic 
arrays,  RF  spatial  arrays,  and  MTI  doppler  filters  [7-11].  Thus,  it  is 
natural  to  consider  their  application  to  ISAR  (Inverse  Synthetic  Array  Radar) 
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imaging,  in  the  hope  that  one  could  obtain  either  higher  resolution  images 
from  the  same  data  samples,  or  equal-quality  images  from  significantly  fewer 
data  samples. 

2.  Doppler  Domain  Processing  Algorithms 

In  this  section,  we  assume  an  input  from  the  radar  receiver  consisting  of 
conventional  range  data  samples  occuring  at  the  usual  Nyquist  sampling  rate, 
f_,  which  is  equal  to  the  IF  receiver  bandwidth.  If  the  radar  utilizes  a 

O 

pulse  compression  waveform,  then  the  sampling  is  performed  after  the  pulse 
compression  filter.  The  spacing  of  the  samples  in  the  time  domain  is  the 
reciprocal  of  f  and  defines  a  "rangebin"  for  the  purposes  of  this 
discussion.  Furthermore,  we  assume  that  a  selectable  "range  window"  of 
adjacent  rangebin  data  samples  is  available  from  each  prf,  and  that  there 
exists  within  this  window  a  prominent  point-target  scatterer  which  serves  as  a 
range  reference  tracking  point,  i.e.,  the  range  window  is  referenced  to  some 
point  on  the  target  and  tracks/moves  with  it.  The  acquisition  and  tracking  of 
this  reference  point  is  cruicial  to  achieving  high  resolution  in  two- 
dimensional  target  imaging,  because  it  is  utilized  for  automatic/adaptive 
self-calibration  of  the  system  [12].  This  reference  point  is  sometimes 
refered  to  as  a  "dominant  scatterer"  on  the  target.  Reference  [12]  is 
recommended  for  a  tutorial  discussion  on  the  need  for  effective  self¬ 
calibration  and  an  adaptive  technique  that  is  capable  of  providing  the 
necessary  compensation. 

In  essence,  the  inverse  synthetic  array  is  "focused"  upon  the  tracked 
dominant  scatterer.  Further  information  on  target  tracking  systems  may  be 
found  in  [13]. 


2 


The  Doppler  domain  is  approached  from  the  viewpoint  of  a  classical  time 
domain  transversal  f ilter[1 3, 1 4] .  as  illustrated  in  Fig.  1.  Here,  we 
visualize  feeding  the  above  rangebin  data  samples  into  a  tapped  delay  line 
wherein  the  inter-tap  delay  is  precisely  equal  to  the  reciprocal  of  the  prf. 
Thus,  the  sequence  of  range  window  data  samples  will  be  synchronized  in  time 

at  every  tap  along  the  delay  line.  With  this  viewpoint  in  mind,  consider 
forming  a  large  digital-data  storage  matrix,  £,  of  dimension  M  rows  by  N 
columns,  wherein  the  columns  contain  the  acjacent  rangebin  data  samples  within 
our  selected  range  window.  Thus,  a  typical  matrix  element,  Dmn,  contains  the 
m1'^  rangebin  sample  in  the  n*'1’1  prf  period.  As  the  data  samples  associated 
with  each  prf  period  corm  in,  two  matrix  processing  steps  occur: 

a)  The  columns  of  _D  are  shifted  to  the  left  by  one,  with  the  oldest 
column  (n=1 )  being  bumped  off  (discarded). 

b)  The  new  incoming  data  samples  from  the  current  prf  are  entered 
into  the  right-hand  column  (n*N). 

Note  that  it  requires  N  prf  periods  to  initially  fill  the  matrix  D_.  After  all 
columns  have  been  filled,  via  the  above  two  steps,  we  can  form  a  signal 
vector,  ^(k),  for  each  rangebin  row, 


E„(k)  =  LE.U), 
— m  ml 


Em2(k)’  Em3(k)’*”’  EmN(k)] 


(1) 


L  L. 

where  (k)  denotes  the  k  n  time  sampling  of  the  signal  vector,  subscript  m 

denotes  the  mth  rangebin,  the  subscripts  1  through  N  denote  the  delay  line 

taps,  and  superscript  t  denotes  the  transpose  of  the  vector  E  (k)  . 

—  m 
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The  time  domain  signal  vector  represented  by  equation  (1)  may  be 
processed  via  many  different  transversal  filter  algorithms  to  obtain  Doppler 


domain  spectra.  Brief  descriptions  of  a  few  such  algorithms  follow. 

Discrete  Fourier  Transform  (DFT) 

The  DFT  is  a  Fourier  representation  of  a  finite-length  sequence  which  is 
itself  a  sequence  rather  than  a  continuous  function  [15].  Not  only  is  the  DFT 
primary  in  conventional  ISAR  image  processing,  but  it  plays  a  central  role  in 
the  implementaton  of  a  variety  of  digital  signal  processing  algorithms, 

including  some  of  the  superresolution  techniques.  One  of  its  implementations 
important  to  our  purposes  is  the  formulation  of  a  frequency  domain  filter 
which  will  permit  sampling  at  a  selected  Doppler  frequency.  Briefly,  we 
multiply  the  signal  vector  of  (1)  with  a  filter  vector  weighted  for  the 
desired  Doppler  frequency, f, 


V(k,f) 

m 


E"(k)-W(f) 
— m  — 


=  z 

n=1 


E  (*)  Wn(f) 

mn  n 


(2) 


where  Vm(k,f)  is  the  kth  time  sample  at  frequency  f  for  the  mth  rangebin. 
_W(f)  is  the  filter  vector  at  frequency  f, 

Wn(f)  is  the  nth  tap  weight  of  W(f). 

Equation  (2)  essentially  defines  a  "beamformer"  as  illustrated  in  Fig.  2.  The 
individual  weights  are  readily  computed  from  the  phase  gradient  along  the 
delay  line  taps  for  frequency  f, 


wn(f) 


expCjuj(N-n)] 


where  w 


radians 


(3) 


(4) 


Uniform  amplitude  weighting  is  employed  in  (3)  and  it  will  result  in  a  (sin 
x)/x  type  of  filter  response.  Tapered  amplitude  weighting  is  an  option,  of 
course. 

Note  in  equation  (4)  that  Doppler  frequency,  f,  is  normalized  to  the  prf  such 
that  our  normalized  angular  frequency,  uj  ,  will  be  in  radians  with  a  range  of 
— it  S  a)  £  it  under  normal  Nyquist  sampling  assumptions  [15]. 

A  second  DFT  implementation  utilized  in  our  processing  is  to  form  a 
multiple  set  of  (sin  x)/x  Doppler  filters  that  satisfy  orthogonal  separation 
criteria.  This  becomes  equivalent  to  a  digital  Fast  Fourier  Transform,  (FFT) 

[15]  and  has  an  exact  correspondence  to  the  Butler  matrix  spatial  beamformer 

[16] .  The  complete  transformation  matrix,  _B,  consists  of  N  column  vectors  of 
N  tap  weights.  A  typical  matrix  element  weight,  b  ^  ,  is  of  the  form, 


(5) 


where  4  =  Doppler  frequency  filter  index 

n  =  tap  weight  index 
N  =  total  number  of  taps 

This  transformation  is  utilized  in  the  "beamspace"  LPF  algorithm. 

A  third  DFT  implementation,  utilized  for  plotting  the  Doppler  spectrum  is 
to  average  the  sampled  power  from  (2)  over  K  time  samples, 


P  (f) 
m 


1 

K 


V*(k,f)*V  (k,f ) 
m  m 


(6) 


where  P  (f)  is  the  Doppler  spectrum  power  at  frequency  f,  and  the  symbol  * 
denotes  complex  conjugate.  Substituting  into  (6)  from  (2)  we  obtain, 
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(8) 


where 


R 

m 


1 

K 


K 

£ 

k=1 


* 

E 

-m 


(k) 


£ 


R  is  the  familiar  sample  covariance  matrix  averaged  over  K  updates  [14]. 
— m 

Note  that  the  total  number  of  prf  implied  by  (8)  will  be  equal  to  the  K 
updates  plus  the  initial  N  prf  required  to  fill  the  delay  line. 


Howells-Applebaum  (H-A)  Algorithm 

This  well-known  adaptive  array  algorithm  is  popular  in  practical 
applications  because  it  does  not  require  data  storage  as  in  the  sample 
covariance  matrix,  nor  does  it  require  a  matrix  inverse  computation.  It 
simply  computes  the  cross-correlations  between  the  array  output  and  each  tap 
signal,  and  integrates  those  updates  into  the  N  weights.  Fig.  2  can  serve  as 
a  schematic  diagram  if  we  consider  the  weights  adaptive.  A  simple  digital 
recursive  version  of  this  algorithm  may  be  written  per  references  [17,18], 


T  W  (k+1 )  =  ( T— 1 )W  (k)  ♦  S  *(f)  -  g  E  *  (k)Y  (k)  (9) 

n  n  n  mn  o 


N 

where  Y  (k)  =  Y'  E  (k)  W  (k)  is  the  array  output  (10) 

o  mn  n 

n=1 


Wn(k+1 )  is  the  updated  n1"*1  tap  adaptive  weight, 
Wn(k)  is  the  current  n  tap  adaptive  weight, 

Emn(k)  is  the  tap  signal  from  equation  (1), 

* 

S  (f)  is  a  deterministic  "steering"  filter  weight, 


g  is  a  gain  parameter, 


T  is  a  weight  integration  parameter. 


The  H-A  algorithm  is  a  steepest-descent ,  feedback  type  of  adaptive  filter  that 
attempts  to  drive  the  output  toward  white  noise  by  moving  filter  zeros/nulls 
onto  spectral  lines  within  the  incoming  signal.  It  requires  only  2N 
multiplies  per  update,  and  yet  is  capable  of  satisfactory  adaptive  filter 
performance  for  many  applications.  However,  it  does  have  certain  performance 
limitations  that  should  be  evaluated  [I1*]. 

Sample  Matrix  Inverse  (SMI)  Algorithm 

The  SMI  algorithm  computes  an  optimum  Weiner  filter  adaptive  weight  from 
the  inverse  of  the  sample  covariance  matrix  [  1  ■h  ] , 

Wq  =  uffV(f)  (11) 

where  Wq  i3  the  optimum  filter  weight  vector, 

* 

S^  (f)  is  a  "steering"  filter  vector  at  frequency  f, 
u  is  a  scalar. 

This  algorithm  is  generally  superior  in  performance  to  the  H-A  algorithm,  but 
it  also  requires  data  storage  in  the  sample  covariance  matrix  and  must  compute 
the  inverse.  Approximately  multiplies  per  update  are  required. 

To  ease  the  computer  storage  requirement  in  accumulating  a  total  of  M  sample 
covariance  matrices,  it  is  expedient  to  store  only  the  upper  triangle  of  the 
Hermitian  matrices,  i.e., 

accumulate  R  = 

— m 
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where  0  denotes  a  lower  triangle  matri::  of  zeros.  When  the  Rm  have  been 

2 

completed,  we  then  compute  and  inject  a  diagonal  "psuedonoise"  term,  6m  , 


where  a2  is  the  average  noise  power  level  per  tap,  and  a  is  the  psuedonoise 
o 

magnitude  factor. 

The  advantages  of  adding  a  small  amount  of  psuedonoise  are  described  in 
reference  C 9 3  - 

The  next  processing  step  then  consists  of  inverting  the  M  modified  sample 
covariance  matrices, 

'  2  -1 

R  =  R  +  B~I  0*0 

— m  I  — m  m- 

i-  J 

where  I  denotes  the  identity  matrix,  and  substituting  into  (11)  to  obtain  the 
optimum  Weiner  filter  weights. 

Maximum  Likelihood  Method  (MLM)  Algorithm) 

The  MLM  algorithm  was  introduced  by  Capon  [6,7],  and  is  defined  as  a 
filter  designed  to  pass  the  power  in  a  narrow  band  about  the  signal  frequency 
of  interest,  and  to  minimize  or  reject  all  other  frequency  components  in  an 
optimal  manner.  Thus,  our  Doppler  spectrum  power  may  be  estimated  by 
substituting  optimum  Weiner  filter  weights  from  (11)  into  (7), 


p  ir) 


R 


-m 


W 


— o 


(15) 
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The  MLM  algorithm  further  imposes  a  unity  gain  constraint  at  each  Doppler 
frequency  such  that, 


3t(f)  W  =  1  06) 

—  — o 

Applying  this  constraint  to  (11)  and  substituting  back  into  (15)  results  in, 

- _J - - -  (17) 

St(f)-R  _1-S*(f) 

—  — m  — 

Therefore,  upon  sweeping  the  deterministic  "steering"  vector,  S^(f),  through 
the  Doppler  frequency  range  of  interest,  the  estimated  power  spectrum  will  be 
plotted  for  the  mth  rangebin.  In  addition  to  estimating  the  spectral  lines, 
the  MLM  algorithm  provides  other  useful  performance  characteristics: 

a)  It  evaluates  the  relative  power  levels  of  estimated  spectral  lines 
because  of  the  gain  constraint  in  equation  (16). 

b)  It  is  one  of  the  few  adaptive  algorithms  that  permits  "copy"  of  a 
desired  narrow  spectral  band  [8],  i.e.,  it  will  pass  the  desired 
narrow  band  while  rejecting  all  other  frequency  components  in  an 
optimal  manner.  This  feature  will  be  used  extensively  in  section  4. 

It's  computation  burden  is  similar  to  that  of  the  SMI  algorithm,  i.e.,  it 
requires  data  storage  in  the  sample  covariance  matrix  and  must  compute  the 
inverse. 

"Beamspace"  LPF  Algorithm 

This  technique  is  patterned  after  a  spatial  beamspace  method  described  by 
Gabriel  [19],  and  is  intended  to  permit  focussing  upon  a  small  sector  of  the 
overall  image.  It  achieves  wide-aperture  high-resolution,  but  yet  requires 


P  (f)  = 
m 
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only  a  relatively  small  number  of  degrees  of  freedom  (DOF)  in  the 
processing.  Fig.  3  illustrates  the  technique  applied  to  a  time-domain  Doppler 
processor,  wherein  a  64-tap  delay  line  receives  the  input  data  samples.  Thus, 
dimension  N=64  in  the  D_  matrix  described  above.  System  DOF  are  reduced  in  two 
steps:  first  a  multiple  sub-band  DFT  connects  into  the  even-numbered  taps 
only  (  a  32-point  DFT);  and  second,  we  select  only  15  of  the  output  sub-bands 
to  cover  our  Doppler  sector  of  interest.  The  resultant  output  Doppler  filter 
passbands  are  shown  plotted  in  Fig.  4.  Fig.  3  also  shows  a  shaped  reference 
filter  connected  into  seventeen  (17)  adjacent  taps,  numbers  48  through  64. 
The  tap  weighting  is  chosen  to  pass  only  the  frequencies  within  the  Doppler 
sector  of  interest,  and  to  reject  all  others.  Fig.  5  illustrates  the  response 
of  our  example  reference  filter.  Thus,  we  end  up  with  just  sixteen  filter 
outputs  to  process,  instead  of  64,  and  can  proceed  with  almost  any  of  the 
current  high-resolution  algorithms,  including  the  algorithms  described  above. 


3.  Simulated  Data  Results 

Rotating  Boom  Model:  Simulated  data  samples  are  generated  from  several 
different  point-target  models.  Fig.  6  illustrates  the  first  of  these,  based 
upon  a  rotating-boom  along  which  are  located  fifteen  (15)  targets.  The  angle 
of  the  boom  with  respect  to  the  direction  of  the  radar  is  defined  as  6 ( t ) , 


9(t) 


(18) 


where  f  is  the  rotation  rate  of  the  boom  in  revolutions  per  minute  (rpm),  and 
t  is  the  time  variable  in  seconds.  The  range  of  the  ith  target,  X^t),  from  a 
zero-range  reference  point  may  be  written, 


10 


(19) 


X. (t)  -  D  -  d.cos(0(t)) 

1  ci 

where  is  the  distance  from  the  center  of  the  boom  to  the  ith  target,  and  Dc 
is  the  distance  from  the  zero-range  reference  point  to  the  center  of  the 
boom.  Range-rate  velocity,  V^(t),  of  the  ith  target  is  given  by  the  first 
derivative  of  X^t), 

dX  <t)  /M 

Vi(t)--dt - 2lt\60/  i3in(  6<  t ) )  i  (20) 

and  we  can  then  readily  compute  the  Doppler  frequency,  f^t),  via  the  RF 
wavelength,  X  , 

fi(t)  •  (!)vi  (t)  (2,) 

The  factor  of  2  in  (21 )  occurs  because  of  the  two-way  path  of  radar  RF  energy 
[13],  Table  1  contains  the  parameter  values  associated  with  the  15  point- 
targets  located  along  the  boom,  together  with  their  Doppler  frequencies 
computed  from  (21).  Five  targets  were  given  a  signal-to-noise  ratio  (SNR)  of 
23  dB,  two  of  20  dB,  one  of  10  dB,  and  the  remainder  0  dB.  Some  targets  were 
grouped  together  to  deliberately  produce  closely-spaced  Doppler  signals  that 
would  require  superresolution  to  separate.  Data  samples  were  accumulated 
from  200  prf  (total)  as  the  boom  swung  through  broadside. 

Complex  data  samples  as  utilized  in  (1)  may  be  computed  for  any  given 
sample  time,  t,  that  is  synchronized  with  the  radar  prf  and  our  selected  range 
window, 
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(22) 


E  ( t )  »  I  A, ( t )  +  n  (t) 
m  1=1  *  ffi 


c  •  c0  *  (577) 


where  tn  is  the  rangebin  index 

Lm  is  the  number  of  targets  within  the  nrn  rangebin 

n  (t)  is  the  receiver  noise  sample  for  the  mth  rangebin 
m 

k  is  the  time  sampling  index 
tQ  is  a  start-time  parameter 

The  A^(t)  are  determined  from  a  rangebin  test  conducted  for  all  pertinent 
targets  at  time  t, 


Vxi(t)  s 


where  A  is  the  rangebin  dimension 

rm  is  the  mth  rangebin  distance 

and  X^(t)  is  defined  in  (19).  Note  that  the  index  l  will  be  associated  with 
our  target  index  i  as  determined  from  (24),  and  that  we  simply  substitute 
those  particular  targets  (if  any)  into  (22)  as 


jr.  (t) 

A£(t)  »  A.e  A 


/x;(t)\ 

,.(t)  . 
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where  the  amplitudes,  Aj,  correspond  to  the  SNR's  given  in  Table  1.  The 
reader  will  note  that  the  target  phases  in  (26)  depend  upon  the  range  in 
wavelengths  at  each  prf  sampling  time  t. 

Fig.  7  shows  the  Doppler  spectrum  estimates  for  nine  of  the  targets  along 
the  boom,  computed  from  the  data  taken  in  the  particular  rangebin  occupied  by 
the  boom  in  its  broadside  position,  i.e.,  this  is  a  Doppler  spectrum  for  a 
single  rangebin  cut.  Utilizing  the  processor  arrangement  illustrated  in  Fig. 
3,  our  delay  line  gives  us  an  overall  Doppler  gradient  window-aperture  of  N=64 
taps,  which  implies  a  sub-band  filter  conventional  resolution  of  360/64  =  5.63 
degrees  in  the  z-plane  (prf-normalized  Doppler  radians  converted  to  degrees) 
Fig.  4  verifies  this  sub-band  filter  resolution.  Therefore,  note  in  fig.  7 

that  we  have  a  target  doublet  and  a  target  triplet  in  which  the  Doppler 
frequencies  are  separated  by  less  than  3  degrees,  or  about  half  of 

conventional  resolution  capability.  Fig.  7(a)  demonstrates  that  the 
conventional  Fourier  transform  is  not  able  to  resolve  the  nine  individual 
point-targets  that  occur  within  the  plot  region.  Yet  those  targets  have  been 
resolved  by  both  of  the  spectral  estimation  algorithms  utilized,  clearly 
demonstrating  superreqolution  for  this  example.  In  Fig.  7(c)  we  note  that  the 
MLM  algorithm  not  only  resolves  the  targets,  but  also  gives  the  correct 
relative  power  level.  Finally,  on  the  right-hand  side  of  the  Fig.  7  plots, 
there  are  four  Doppler  signals  of  significantly  different  power  levels  and 
spaced  one  beamwidth  apart.  This  type  of  situation  defies  resolution  by 
conventional  filters  of  64  taps  because  of  their  beamwidths  and/or  sidelobes, 
and  yet  the  two  superresolution  algorithms  have  readily  resolved  the  four 
signals.  A  conventional  FFT  processor  would  require  an  aperture  several  times 
larger  than  our  64  taps  in  order  to  resolve  these  nine  individual  point- 
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targets,  and,  therefore,  may  actually  require  more  prf  data  samples  to  achieve 
the  same  resolution.  Note  in  this  example  that  N=64,  K*140  updates  per  (8), 
and  there  was  a  total  of  200  prf  data  samples  available. 

A  point  of  clarif ication  is  in  order  concerning  the  appropriate  measure  of 
resolution  for  comparison.  It  is  well  known  that  temporal  frequency 
resolution  varies  inversely  wi' h  the  processing  time  duration  [5].  For  the 
example  given  in  TABLE  1,  we  have  a  prf  of  418  and  200  samples  processed, 
resulting  in  a  processing  time  duration  of  0.48  seconds  or  a  temporal 
frequency  resolution  of  2.1  Hz.  However,  this  ultimate  resolution  assumes 
uniform  weighting  of  the  data  samples  plus  the  utilization  of  all  200  samples 
to  get  a  single  output  pulse  in  adjacent  filter  bins  -  assumptions  that  are 
often  not  applicable  in  practical  systems.  Doppler  filters  often  have  tapered 
weighting  to  reduce  sidelobes  and,  in  addition,  utilize  some  form  of  output 
pulse  integration.  Thus.,  a  reasonable  measure  for  resolution  comparison 
should'  be  based  upon  the  particular  system  under  consideration. 

Targets  Entering/Leaving  a  Rangebin:  Before  we  begin  to  demonstrate  three- 
dimensional  (3D)  imaging  of  rotating  objects,  it  is  necessary  to  address  the 
deleterious  effects  associated  with  targets  that  are  either  entering  or 
leaving  a  rangebin.  Whenever  the  returns  from  a  target  fail  to  fill  our 
tapped  delay  line  in  a  given  rangebin,  then  we  no  longer  have  a  complete 
gradient  aperture  and  the  target  energy  will  not  focus  into  a  narrow  Doppler 
banci.  Instead,  the  energy  tends  to  spread  out  in  accordance  with  its  partial 
aperture,  resulting  in  a  smeared/ faulty  Doppler  spectrum  estimate. 

One  technique  for  diminishing  such  faulty  spectrum  effects  is  to  monitor 
the  tap  power  in  both  halves  of  the  delay  line  and  form  two  sums  that  would  be 
sensitive  to  a  strong  target  entry-exit, 
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-  \  w*! 

n- 1  1 


left  half  sum. 


S  (k)  =  Z  E  (k) 
mr  ,  mn 

n«L 


,  right  half  sum 


where  L 


and  Emn(k)  is  defined  in  (1), 


These  sums  must  usually  be  augmented  with  a  bias  term  to  prevent  excessive 
noise  triggering. 


si  ■  s„t(k) 


s«r‘k)  11  «o 


where  aQ  is  the  average  noise  power  level  per  tap.  We  then  form  a  power- 
sensitive  threshold  level,  Q,  which  may  be  of  the  form, 


Q  =  1.81  + 


If  the  ratio  of 


ft) 


V  s2 


or  its  reciprocal,  exceeds  the  value  of  Q,  then  that 


particular  rangebin  signal  vector  is  skipped,  i.e.,  it  is  not  utilized  in  the 
algorithm  update  processing. 

A  better  technique  that  is  applicable  to  targets  of  low  SNR  is  to  apply 
an  FFT  to  both  halves  of  the  delay  line,  such  that  we  now  divide  the  rangebin 
signal  vector  energy  into  L  Doppler  filter  bins  as  indicated  in  (2)  through 
(5).  Note  that  each  filter  bin  has  a  power  gain  equal  to  L,  thus  enhancing 
targets  of  low  SNR.  If  a  target  is  entering  a  rangebin,  its  Doppler  energy 
will  be  strong  in  its  right-half  Doppler  bin,  but  non-existent  in  the  same 
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Doppler  bin  on  the  left-hand  side;  and  vice-versa  if  a  target  is  leaving  a 
rangebin.  Therefore,  the  ratio  of  the  power  in  each  Doppler  bin  must  be 
compared  for  the  two  halves  in  a  manner  very  similar  to  (29)  through  (31),  and 
particular  rangebin  signal  vectors  are  skipped  if  the  threshold  is  exceeded. 
The  performance  of  this  FFT  technique  is  much  superior  to  the  simpler 
technique  described  in  (27)  through  (31 ),  but  it  also  requires  considerably 
more  computation  in  doing  the  FFT's  and  then  testing  each  Doppler  filter  bin. 
Merry-go-round  Model:  The  second  simulation  model  is  based  upon  a  "merry-go- 
round"  of  24  point-targets,  as  illustrated  in  Fig.  8.  The  targets  are  equally 
spaced  on  a  circle  of  35  meters  radius  and  the  entire  assembly  is  rotated  at  a 
rate  of  0.63  rpm.  Target  reflection  coefficients  are  chosen  such  that  the  odd 
numbers  have  a  receiver  output  signal-to-noise  ratio,  (SNR)  of  5  dB,  and  the 
even  numbers  have  an  SNR  of  25  dB.  Thus,  they  alternate  in  strength  around 
the  circle.  In  addition,  a  stationary  target  of  25  dB  SNR  is  placed  at  the 
center  of  the  circle  to  serve  as  a  zero-Doppler  reference. 

From  Fig.  8,  we  can  compute  the  angle  position  of  the  ith  target,  ei  (t), 
at  time  t 


e.(t) 


(32) 


where  I  is  the  total  number  of  equi-spaced  targets  and  fp  is  the  rotation  rate 
in  rpm.  The  range  of  the  it*1  target,  X^(t),  from  a  zero  range  reference  point 
may  be  written, 

Xi(t)  =  Dc  -  docos(8i(t))  (33) 
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where  Dc  is  the  distance  from  zero-range  reference  point  to  the  center  of  the 
circle,  and  dn  is  the  radius  of  the  circle.  Range-rate  velocity,  V^t),  and 
Doppler  frequency,  f^(t),  then  follow  as  in  (20)  and  (21).  Complex  data 
samples  may  be  generated  from  this  model  in  the  same  manner  as  described  for 
the  rotating  boom  model. 

Fig.  9  shows  a  typical  three-dimensional  range-Doppler  radar  image  for 
this  particular  rotating  object.  The  upper  portion  of  Fig.  9  illustrates  the 
"image  truth",  i.e.,  it  plots  the  point-target  SNR's  at  their  exact  range- 
Doppler  coordinate  locations  on  the  "merry-go-round" .  The  lower  portion  of 
Fig.  9  illustrates  a  3D  image  estimate  of  the  rotating  object,  utilizing  the 
simple  H-A  algorithm  described  in  the  previous  section.  Parameter  values 
include  M=90  rangebins,  N=32  delay  line  taps,  and  K  =  140  prf  processed. 
Comparing  the  estimate  against  the  "image-truth" ,  we  note  that  the  H-A 
algorithm  has  correctly  located  all  24  point-targets  around  the  circle.  On 
the  negative  side,  this  simple  algorithm  does  not  give  us  the  relative  power 

levels  of  the  target,  and,  in  addition,  it  has  a  rather  "noisy  floor"  that 
will  require  thresholding  to  enhance  the  image. 

Fig.  10  illustrates  a  second  3D  image  estimate  of  this  rotating  object, 
utilizing  the  more  desirable  MLM  algorithm  described  in  the  previous 
section.  Recall  that  the  MLM  algorithm  has  a  much  greater  computation  burden 
than  the  H-A  algorithm  but,  in  return  it  also  estimates  the  relative  power 
levels  of  the  targets  in  addition  to  their  locations.  Fig.  10(a)  was 
processed  without  any  regard  for  target  entry/exit  effects,  and  it  contains 
several  examples  of  the  smeared/faulty  Doppler  spectrum  problems  referred  to 
above.  In  contrast,  Fig.  10(b)  was  processed  utilizing  the  FFT  entry/exit 
routine,  and  it  will  be  noted  that  all  of  the  problem  target  spectra  were 
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eliminated.  Another  technique  utilized  in  Fig.  10  to  enhance  the  image  was 
the  use  of  a  2  dB  output  threshold  level,  which  tends  to  diminish  the  "noisy 
floor"  effect  (compare  against  Fig.  9(b)). 

A.  Range  Domain  Superresolution 

The  problem  of  resolving  an  unknown  number  of  closely  spaced, 
overlapping,  and  noisy  echoes  of  a  signal  with  apriori  known  shape  is  common 
to  several  application  areas  such  as  radar,  sonar,  geological  acoustic 
sounding,  ultrasound-based  nondestructive  testing,  and  medical  imaging 
procedures.  Several  approaches  for  solving  this  problem  have  been  described 
xn  the  literature.  These  include  detection/ deconvolution  schemes,  inverse 
filtering,  least  squares,  and  maximum  likelihood  methods  [20-25].  In  the 
radar  area,  we  are  familiar  with  various  pulse  compression  schemes  [133  that 
have  been  applied  to  this  problem.  Recently,  Bruckstein  et  al  [26]  have 
described  a  new  method  for  achieving  high  resolution  in  such  cases,  based  upon 
an  eigenstructure  technique  which  exploits  the  structure  of  the  received 
signal  covariance  matrix.  This  new  method  was  instrumental  in  prompting  the 
current  investigation.  Also,  high  resolution  techniques  have  been  applied  to 
the  wideband  linear  FM  ranging  method  [5,27,28]. 

Consider  the  familiar  rectangular  pulse  with  linear  FM  chirp  shown  in 
Fig.  11.  This  type  of  pulse  [133  is  utilized  as  an  example  of  a  radar  signal 
of  known  shape,  S(t),  throughout  the  remainder  of  our  discussion.  S(t)  may  be 
expressed  mathematically  in  the  form, 

S(t)  =  rect^|  ^  exp^^-jr^  \  -  |  <  t  <  |  ,  (3*0 
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where  T  is  the  pulse  duration  (width),  B  is  the  modulation  bandwidth,  and  t  is 
the  time  variable.  The  exponent  term  gives  us  the  phase  modulation  during  the 
pulse,  'f(t), 

Hf(t)  =  ’  -  \  <  t  <  |  (35) 
and  the  first  derivative  of  't'(t)  gives  the  frequency  modulation  (chirp) 


(36) 


Characteristics  of  interest  include  the  time- bandwidth  product,  BT,  the 
approximate  time  resolution  given  by  the  reciprocal,  ^  ,  and  the  pulse 

compression  ratio,  D, 


D  =/T\  -  BT  (37) 

w 

We  are  also  interested  in  delayed  and  scaled  versions  of  our  known  waveform 
that  would  correspond  to  the  range  locations  of  scattering  objects/ point 
targets, 


a.S(t-T. ) 

■k  1 


(T.-  f)  <  t  «Tt.  \) 


where  -  (——)  ,  the  round-trip  time  delay  associated  with  the  ith  point- 
target  at  range  rif  c  is  the  velocity  of  light,  and  is  the  scattering 
coefficient. 
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The  assumptions  discussed  in  the  beginning  paragraph  of  Section  2  also 

apply  here,  except  for  one  important,  fundamental  difference.  For  our  range 

domain  processing,  the  sampling  must  be  performed  prior  to  the  pulse 

compression  filter.  The  reason  for  this  requirement  is  that  the  pulse 

compression  filter  will  be  utilizing  adaptive  weighting,  rather  than 

conventional  deterministic  matched  filter  weighting.  One  of  the  implications 

of  this  fundamental  difference  is  that  we  must  first  define  an  uncompressed, 

digital-data  signal  vector,  U  ,  within  our  range  window, 

k 


W  *  W’  W’*"’  uk(V 


) 


(39) 


where  the  subscript  k  denotes  the  kLn  prf  period  ,  and  subscripts  1  through  Q 
denote  the  q  n  data  sample.  This  uncompressed  signal  vector  may  be  processed 
via  various  pulse  compression  algorithms  to  obtain  desired  range-domain 
response.  Brief  descriptions  of  a  few  such  algorithms  follow. 

Matched  Filter  DFT 

First,  let  us  compute  the  conventional  matched  filter  pulse  compression 
response  wherein  we  require  matched  filter  weights,  W(t)  at  time  t=*  T 


corresponding  to  a  given  range.  Briefly,  we  desire  the  complex  conjugate  of 
(3*0,  but  including  a  linear  phase  term  to  account  for  the  time  delay  phase 
shift, 

W(t-T)  «  recti ^—jexp  (-j'i'(t-T))  (40) 


and  ¥(t-T)  =  ^pj(t-T)2  -  u»qT  (41) 

subject  to  (T  -  |  )  <  t  <  (T  +  |  ) 
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where  ojqT  is  the  linear  phase  term  to  account  for  the  time  delay  phase 

shift.  In  practice,  there  is  a  total  linear  phase  term  of  u>  (t-T),  but  we 

delete  the  reference  RF  carrier  phase,  u  t  ,  for  the  sake  of  simplicity  in  the 
expressions.  Since  we  are  dealing  with  a  precise  sampled-data  system,  (40) 
and  (41)  must  be  written  in  terms  of  integer  multiples  of  the  sampling 
time,  At  ,  which  is  hereafter  assumed  to  be  normalized  to  unity,  i.e., 
At  =  — g  -  1 .  This  implies  that  our  bandwidth,  B,  is  also  normalized  to 

unity.  Further,  we  may  account  for  the  carrier  time  delay  phasing  via  the 
number  of  RF  wavelengths  per  At  interval. 

Note  that  our  signal  vector  in  (39)  corresponds  to  this  model,  and  we  have 
these  precise  sampling  times  denoted  as  t^,  where  the  range  window  extends 
from  1<q<Q.  The  waveform  of  Fig.  11  has  a  width  smaller  than  Q  and,  from 
(37)  we  find  that  the  number  of  samplings  ,N,  within  the  waveform  is  equal  to 

the  pulse  compression  ratio,  i.e.,  define  N  *  INT(D)  where  INT(  )  means  the 

integer  value  of  the  argument.  Applying  the  waveform  time  restriction,  it 
follows  that  the  values  of  within  the  waveform  for  delay  T  must  be, 


or 


(T-  h  <  t  < (T+~ ) 
2  q  2 


-  \  <(yT)<I 


(42) 


If  T  is  expressed  in  terms  of  sampling  units,  then  we  may  define  an  integer 
index  1  ,  such  that 

l  -  INT(T+.5)  (43) 
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and  this  integer  serves  as  an  index  for  individual  rangebins  within  the 
overall  range  window.  Also,  it  permits  us  to  define  an  index  n  that  applies 
only  within  the  waveform, 


define  t 

n 


rj(2n-N-1)  ♦  Z 


for  N  odd, 


(44) 


such  that 


(tn  '  T) 


(#5) 


Therefore  (t  -  T)  »  ~  (2n-N-1)  -  (T-l)  (46) 

where  1  £  n  2  N 

From  (42)  it  is  evident  that  the  t^  may  be  associated  with  N  values  of  t^  via 
a  simple  index  shift, 

q  *  n  +  l  {Hi) 


Substituting  (46)  into  (40)  and  (41)  gives  us  our  matched  filter 
weights,  W  (t  -T)  .and  we  can  apply  them  to  the  signal  vector  Uj,  to  obtain  a 

conventional  OFT  compressed  pulse  output, 

N  , 

V  (T)  -  Vujt  )  W  (— (2n-N-1  )  -  (T-U)  (48) 

k  ,  k  q  n  2 

n-1  ^ 

where  index  q  is  given  by  (47)  above.  The  variable,  T,  becomes  our  parameter 

2r 

for  "sweeping"  or  "stepping"  across  the  range  window.  Recall  that  T-(— ). 

Equation  (48)  illustrates  mathematically  how  the  waveform  operates  upon 
the  signal  vector  to  generate  a  DFT  compressed  pulse  output  for  a  target 
located  within  a  particular  rangebin.  Note  that  the  variable,  (T  -l)  , 
functions  as  a  vernier  shift  for  each  of  the  N  sampling  times  within  the 
waveform. 
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MLM  Algorithm 


The  description  of  this  algorithm  in  Section  2  applies  equally  to  the  range 
domain,  and  it  is  an  important  algorithm  because  it  allows  us  the  "copy" 
privilege.  That  is,  it  can  pass  the  signal  in  a  narrow  band  around  the  target 
range  of  interest,  and  minimize  or  reject  all  other  target  signals  in  an 
optimal  manner. 

To  utilize  the  algorithm,  we  must  form  the  familiar  sample  covariance 
matrix,  R  ,in  the  range  domain.  R  may  be  averaged  over  K  updates  in  a 
manner  similar  to  equation  (8)  in  the  Doppler  domain, 


R 


l 

K 


K 


E 

k*1 


# 


(49) 


However,  an  important  difference  is  that  we  do  not  require  filling  all 
elements  within  the  matrix  dimensions  of  QxQ,  since  the  compression  filter 
weights  only  have  a  dimension  equal  to  N.  Thus,  the  indices  of  our  matrix 
elements,  Rmn,  in  the  stored  upper  triangle  (12)  may  be  limited  to  (n-m+1)£N 
where  m  is  the  row  index,  and  n  is  the  column  index. 

Because  the  dimension  N  is  considerably  smaller  than  Q,  it  is 
advantageous  to  extract  from  R  a  series  of  smaller  matrices  denoted  as  R^  , 
where  their  dimensions  are  exactly  NxN  and  the  subscript,  i,  is  related  to  the 
row/ column  subscripts  as  follows, 

(i+1 )  S  m  S  (i+N)  (50) 

i=0, 1  ,2,3 . (Q-N) 

(i-1 )  S  n  S  (i+N) 
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This  results  in  a  total  of  (Q-N+1 )  smaller  matrices  which  are  to  be  inverted 


and  stored  for  use  in  processing. 

To  compute  the  MLM  optimum  adaptive  weights  for  "copy"  purposes,  we  first 
substitute  (40)  into  (11)  to  get  the  Weiner  weights, 

W  (t  -  T)  =  u  R"1  W(t  -T)  -  5  <  (t  -  T)  <  |  (51) 
—  o  -l  -  2  2 

where  index  i  =  l  per  equation  (43). 

It  should  be  noted  that  we  are  using  the  conventional  matched  filter  vector  as 
our  "steering"  vector.  The  scalar  y  is  next  evaluated  as  the  MLM  power 
output  per  (17),  and  substituted  back  into  (51)  to  give  us  gain-normalized  MLM 
optimum  "copy"  weights.  The  reader  should  recognize  that  (51)  results  in 
adaptive  pulse  compression  filter  weights  which  tend  to  null  out  any  target 
returns  that  are  not  at  the  range  represented  by  T  . 

A  sjttle  requirement  embodied  within  (51)  is  the  small  matrix  index,  i. 
As  we  "sweep"  or  "step"  our  steering  vector  via  the  value  of  T,  then  the 
correct  ith  matrix  inverse  must  be  utilized.  This  index  value  follows  from 

(43)  and  (47).  Each  rangebin  requires  its  associated  small  matrix,  i.e.,  if 
we  are  sweeping  across  a  range  window  consisting  of  30  rangebins,  then  we  need 
30  small  matrices.  A  given  small  matrix  permits  optimum  weights  to  be 
computed  anywhere  within  the  extent  of  its  associated  rangebin,  i.e.,  within  ± 
1/2. 

The  MLM  adaptive  "copy"  weights  are  utilized  in  exactly  the  same  manner 
as  conventional  matched  filter  weights  in  (48).  However,  the  vernier  shift 
permitted  by  (T-H)  is  now  of  considerable  interest  because  of  our  intention  to 
achieve  range  resolution  finer  than  a  single  sampling  period  (or  rangebin). 
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SLC  Algorithm 


This  algorithm  is  derived  from  its  well-known  counterpart  in  the  spatial 
domain,  the  "sidelobe  canceller"  (SLC)  [ 1 ^ ,  17]  which,  in  turn,  is  intimately 
linked  to  the  familiar  linear  prediction  filter  of  spectral  analysis  [6].  The 
algorithm  is  widely  used  because  it  is  relatively  simple  to  implement  in  most 
applications,  and  yet  exhibits  excellent  adaptive  filter  performance 
characteristics.  The  unique  feature  of  this  algorithm  is  its  simple 
"steering"  weight  vector,  which  consists  of  zeros  except  for  one  element.  For 
example,  in  (51)  our  "steering"  vector  becomes, 

Wfc(t  -T)  =  [1,0,0 . 0.0  ,]  (52) 

and  it  is  evident  that  the  optimum  Weiner  weight  will  then  consist  of  the 
first  column  of  the  inverse  of  the  small  sample  covariance  matrix.  To 
identify  its  special  nature,  the  resulting  optimum  weight  is  referred  to  as 
the  Weiner  SLC  adaptive  weight,  and  it  is  standard  practice  to  compute  a  value 
for  u  such  that  the  first  element  of  this  SLC  weight  is  always  normalized  to 
unity. 

A  particularly  attractive  version  of  the  SLC  algorithm  results  when  it  is 
based  upon  the  H-A  algorithm  discussed  in  Section  2,  because  then  the 
computation  burden  is  minimal.  As  explained  therein,  sample  covariance 
matrices  are  not  required,  since  updates  are  integrated  directly  into  the 
adaptive  weights  themselves. 

Coherent  Target  Scatterers 

When  target  objects  contain  scatterers  that  have  no  relative  motion  when 
viewed  from  the  direction  of  the  radar,  then  we  have  a  special  situation 
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defined  as  "coherent  scatterers." 


This  phenomenon  has  an  equivalence  in  the 


spatial  domain  and  has  been  addressed  in  the  literature  by  a  number  of  authors 
[6,11,29].  It  can  result  in  severely  degraded  resolution  performance  from 
superresolution  algorithms,  and  usually  requires  additional  processing 
techniques  to  overcome. 

To  begin  with,  the  reader  is  reminded  that  the  sample  target  data  from  a 
single  radar  waveform  transmission  (prf)  is  inherently  coherent,  in  a  manner 
analogous  to  taking  a  single  data  "snapshot"  across  a  spacial  domain  array  of 
elements.  For  this  singular  situation  the  FFT  matched  filter  (beamformer)  is 
perfectly  operable  at  full  aperture  resolution,  whereas  superresolution 
techniques  are  not.  This  is  one  of  the  reasons  why  superresolution  techniques 
must  be  considered  as  complementary/ junior  to  the  FFT.  Superresolution 
techniques  do  not  become  competitive  until  after  a  reasonable  number  of  prf 
data  samplings  are  available  for  processing. 

The  crucial  r-  "irement  is  to  have  enough  samplings  (observation  time)  to 
permit  decorrela  .on  of  multiple  point-targets  within  a  rangebin  (analogous  to 
multiple  sources  within  a  beamwidth  in  the  spatial  domain).  In  essence,  we 
seek  stationarity  of  the  scatterer  phase  data,  as  if  they  were  sources  of 
random  phase  signals.  There  are  several  ways  in  which  the  necessary 
decorrelation  might  be  achieved,  or  at  least  approximated: 

a.  A  slowly  rotating  object  wherein  the  scatterers  have  sufficient 
doppler  differences.  This  is  a  classic  ISAR  situation  [1-5]. 

b.  An  object  in  straight-line  motion  (non-radial)  wherein  an  equivalent 
rotation  occurs  over  the  observation  time  [5,12]. 

c.  A  stepped  linear  FM  carrier  shift  (small  increment  per  prf)  which 
produces  sufficient  phase  changes  between  scatterers  spaced  closely  in 
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range.  This  borrows  from  the  classic  linear  FM  ranging  method  [5,13]. 

d.  Random  frequency  hopping  of  the  carrier  per  prf  to  produce  sufficient 
phase  changes  for  small  range  separations. 

e.  Combinations  of  the  above. 

It  is  sometimes  helpful  to  view  the  above  decorrelation  methods  in  terms  of 
the  variance  of  the  first  derivative  of  the  ratio  (-)  ,  where  x  is  the  range 

A 

separation  between  two  closely-spaced  scatterers. 

Variance  of  — 1 - £  1  (53) 

dt 

1  H  Y 

Doppler  Difference  =  (— r— )  ^  Hz  (54) 

X  dt 

Carrier  Shift  =  x  d ■  ^-7—  Hz  (55) 

dt 

There  are  a  number  of  situation  conditions  and  suttle  limitations  involved  in 
the  above  methods,  which  are  not  addressed  in  this  initial  report  because  of 
space  considerations.  Further  research  investigations  are  needed  to  clarify 
some  of  these.  The  simulation  examples  contained  in  this  report  are  based 
primarily  upon  method  (a),  the  slowly  rotating  object. 

Simulated  Data  Results 

Simulated  data  samples  have  been  generated  from  a  rotating-boom  model 
similar  to  Fig.  6,  along  which  are  located  three  targets.  TABLE  2  contains 
the  parameter  values,  together  with  their  Doppler  frequencies  computed  from 
(21).  The  parameter  values  were  deliberately  chosen  to  produce  close  target 
spacings  that  would  require  superresolution  in  both  the  range  domain  and  the 
Doppler  domain.  The  pulse  compression  ratio  of  our  chirp  waveform  is  D-1 3 - 
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Uncompressed  digital  data  samples  are  entered  into  our  signal  vector, 

per  (39),  and  are  utilized  in  accumulating  a  sample  covariance  matrix  over 
K*264  prf  updates  per  (49).  Fig.  12a  illustrates  typical  conventional  DFT 
matched  filter  responses  plotted  during  this  accumulation  time,  and  represents 
examples  of  the  V^(T)  output  sequences  per  (48).  The  reason  for  the  large 
spread  in  the  power  plot  values  is  the  continuously  changing  phase 
relationships  between  the  three  targets  due  to  the  different  Dopplers.  Fig. 
12b  shows  the  averaged  matched  filter  response  (7)  over  all  K  prf  updates. 
Fig.  12  clearly  demonstrates  that  our  three  targets  are  too  close  together  to 
be  resolved  by  conventional  pulse  compression.  The  somewhat  "jerky"  behavior 
in  the  range  "sidelobes"  is  due  to  the  waveform  entry/exit  effects  as  it 
sweeps  across. 

After  R  is  accumulated,  we  extract  the  series  of  smaller 
matrices,  per  (50)  and  compute  their  inverses  for  use  in  (51).  Figure  13 
then  illustrates  the  target  estimate  plots  obtained  from  the  MLM  algorithm 
(17),  and  the  Weiner  SLC  algorithm  weight  per  (51)  &  (52).  Note  that  both  of 
these  superresolution  algorithms  clearly  resolve  the  three  targets  and 
accurately  locate  their  range  positions.  In  addition,  the  MLM  algorithm 
correctly  estimates  the  relative  power  levels  of  the  targets. 

The  next  step  is  to  take  advantage  of  the  "copy"  feature  referred  to  in 
the  description  of  the  MLM  algorithm.  This  feature  is  so  important  to  our 
processing  that  it  is  re-iterated  here  again;  i.e.,the  MLM  optimum  "copy" 
weights  in  (51)  permit  us  to  pass  the  signal  in  a  narrow  region  around  the 
target  range  of  interest,  and  minimize  or  reject  all  other  target  signals  in 
an  optimal  manner.  This  feature  allows  us  to  obtain  adaptive  compressed-pulse 
output  signals  at  vernier  steps  in  range.  For  the  current  simulation  example 
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described  in  TABLE  2,  each  sampling  period  (or  rangebin)  is  divided  into  ten 
subdivisions,  and  an  output  data  sample  is  computed  for  each  subdivision  via 
(51)  and  (48).  Fig.  14  shows  four  typical  "copy"  weight  responses  in  the  range 
domain  for  various  "steering"  range  selections.  Starting  at  the  upper  left, 
the  range  mainlobe  is  steered  to  range  12.0  whereby  the  three  targets  are  all 
located  in  the  sidelobe  region  and  readily  nulled  out.  In  the  upper  right 
figure,  the  range  mainlobe  is  steered  to  range  14.0  whereby  the  target  located 
at  14.3  is  encroaching  into  the  mainlobe  and  results  in  high  adaptive 
sidelobes  (compare  with  upper  left  pattern).  However,  note  that  all  three 
targets  are  still  nulled  out.  The  two  lower  figures  illustrate  the  mainlobe 
steered  directly  into  target  locations  at  ranges  14.3  and  15.0.  Note  that  the 
unity  gain  constraint  (16)  results  in  passing  the  target  signal  at  the  steered 
range,  while  simultaneously  nulling  out  the  other  two,  i.e.,  a  clear 
demonstration  of  the  "copy"  feature.  Fig.  15  illustrates  the  results  obtained 
from  applying  this  feature  to  the  rangebin  region  extending  from  range  12.0  to 
18.0.  Compare  against  Fig.  12.  Since  we  are  increasing  the  number  of  output 
data  samples  by  a  factor  of  ten,  it  becomes  expedient  to  apply  this  technique 
over  a  smaller  field  of  view.  Thus,  our  range  field  is  restricted  to  six 
rangebins  in  this  case,  from  range  12.0  to  18.0,  but  that  generates  60  vernier 
data  samples  per  prf.  The  60  data  samples  per  prf  are  then  processed  in  the 
Doppler  domain  in  exactly  the  same  manner  as  described  in  Section  2. 

5.  Combined  Range/Doppier  3D  Superresolution  Imaging 

Three  different  processing  techniques  are  utilized  to  demonstrate 
combined  Range/Doppler  3D  superresolution  radar  imaging: 
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a.  Conventional  Fourier  Transform 


b.  MLM/LPF  Algorithm 

c.  Digital  image  from  optimum  weight  roots. 

Simulated  radar  data  samples  were  generated  from  the  rotating-boom  model 
described  in  Section  4  above,  wherein  we  have  three  point-targets  of  10dB  SNR 
located  along  the  boom.  TABLE  2  contains  the  various  parameter  values  which 
were  deliberately  chosen  to  produce  close  target  spacings  of  0.7  units  in  both 
the  range  and  the  Doppler  domain,  i.e.,  the  targets  are  separated  by  only  0.7 
of  a  conventional  rangebin  in  the  range  domain,  and  0.7  of  a  conventional 
"Doppler  bin"  in  the  Doppler  domain.  Note  the  slow  boom  rotation  rate  of  only 
0.1  rpm. 

Conventional  Fourier  Transform 

Figure  16(c)  shows  the  conventional  Fourier  transform  image  processed  from  the 
simulated  radar  data  samples.  The  Fourier  transform  serves  as  a  baseline 
reference  to  illustrate  what  conventional  processing  would  result  in.  Uniform 
illumination  matched- filter  weighting  was  utilized  in  both  domains,  resulting 
in  maximum  conventional  target  resolution  but  with  the  associated  high 
sidelobes.  An  image  threshold  level  of  15  dB  below  peak  power  level  was 
imposed  in  order  to  eliminate  the  distracting  effects  of  the  high  sidelobes, 
and  it  will  be  noted  that  only  the  first  sidelobes  (13  dB  down)  are  visible  in 
the  image.  At  maximum  conventional  resolution,  we  note  that  the  3D  image  just 
barely  manages  to  resolve  the  three  point-targets  because  their  separation 
along  the  range/ doppler  diagonal  is  unity,  i.e.,  the  square-root  of  the  sum  of 
the  squares  of  their  0.7  spacings.  Recall  from  Fig.  12b  that  resolution  in 
only  the  range  domain  was  not  possible.  Furthermore,  if  we  were  to  utilize  a 
tapered  matched-filter  weighting  in  order  to  reduce  sidelobe  levels,  then  the 
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3D  image  would  no  longer  resolve  the  targets  because  of  the  increase  in  the 
"beamwidths" . 

MLM/LPF  Algorithm 

Fig.  16(b)  shows  the  image  obtained  when  the  same  simulated  radar  samples 
are  processed  via  the  MLM/LPF  Algorithm.  Recall  that  this  algorithm  was 
described  in  Section  2  for  the  Doppler  domain  (see  Fig.  7c),  and  in  Section  4 
for  the  range  domain  (see  Fig.  136).  The  "copy"  feature  of  the  MLM  algorithm 
was  utilized  for  each  vernier  range  subdivision  within  the  small  field  of  view 
represented  by  the  rangebin  scale  in  Fig.  16,  i.e.,  90  vernier  data  samples 
per  prf.  The  MLM  algorithm  has  such  well-behaved  filter  residue 
characteristics  (an  adaptive  equivalent  of  "siSelobes")  that  the  image 
threshold  level  can  be  set  just  slightly  above  the  reciever  noise  level,  about 
2  dB  for-  Fig.  16(b),  thus  eliminating  the  "noise  floor"  from  the  image.  Note 
that  the  three  targets  are  clearly  resolved  in  this  image,  their  locations  are 
correct,  and  their  relative  power  levels  are  accurately  represented.  It  is 
obvious  that  the  three  point-targets  could  be  located  much  closer  together  and 
still  be  resolved.  Although  the  MLM  algorithm  worked  well  in  this  particular 
example,  it  has  a  few  drawbacks  that  must  be  considered: 

a.  The  computation/ processing  burden  is  high. 

b.  It  has  the  poorest  resolution  capability  of  all 
superresol uton  techniques. 

c.  There  is  a  wide  spectral  spread  in  the  Doppler  energy  for 
vernier  range  subdivisions  close  to  a  target  location.  This 
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Doppler  energy  spread  is  caused  by  the  high  sidelobes  of  the 
"copy"  weights  when  our  mainlobe  is  steered  close  to  target 
locations,  i.e.,  the  high  sidelobes  amplify  the  power  level  of 
noise  and  signal  residue  across  the  Doppler  band.  It  produces  a 
"corrugation  effect"  in  the  images,  which  is  evident  in  Fig.  16. 
Digital  Image  from  Optimum  Weight  Roots 

There  are  several  different  methods  available  for  addressing  the 
drawbacks  of  the  MLM  algorithm,  and  the  one  chosen  for  our  next  image  display 
is  described  as  a  "digital"  3D  plot  technique  based  upon  optimum  weight 
roots.  Fig.  16(a)  ‘illustrates  the  3D  image  formed  when  the  same  simulated 
radar  data  samples  are  processed  via  this  technique.  Note  that  the  image  now 
exhibits  excellent  fidelity  with  the  model  truth,  and  there  are  no  distracting 
artifacts  except  for  some  minor  "noise  floor"  pips. 

The  processing  begins  in  the  range  domain  by  accumulating  the  sample 
covariance  matrix  and  extracting  the  series  of  smaller  matrices  denoted  as 
FL  as  described  in  Section  4.  It  is  only  necessary  to  extract  the  nine 

smaller  matrices  associated  with  the  rangebins  to  be  displayed.  These  nine 
matrices  are  inverted,  stored,  and  utilized  to  compute  the  nine  optimum  weight 
vectors  needed,  i.e.,  we  compute  one  adaptive  filter  weight  vector  per 
rangebin.  Fig.  13a  is  a  "range  spectral  estimate"  obtained  from  those  weight 
vectors,  with  the  nine  rangebins  extending  from  range  10.5  to  19.5  as  plotted 
there.  The  sharp  peaks  correspond  to  roots  of  the  adaptive  filter  weights 
which  represent  the  target  range  locations,  so  we  solve  for  those  range 
locations  by  utilizing  a  peakfinder  routine. 

The  peakfinder  routine  readily  detects  the  three  target  locations  within 
its  field  of  view  extending  from  range  10.5  to  19.5.  For  our  current 
simulation  example,  those  three  locations  were  estimated  as  14.3,  15.0,  and 
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15.7.  Knowing  those  locations  gives  us  the  opportunity  for  dramatic 
reductions  in  computer  processing  burden,  because  we  now  require  the  "copy" 
feature  of  the  MLM  algorithm  only  at  the  three  estimated  target  locations, 
instead  of  the  90  vernier  range  divisions  that  would  otherwise  be  computed  to 
cover  the  entire  field  of  view. 

Next,  the  "copy"  vernier  data  samples  for  the  target  locations  are 
processed  in  the  Doppler  domain  via  the  "Beamspace"  LPF  algorithm  of  Section 
2,  from  which  we  compute  associated  adaptive  weight  vectors  for  the  tapped 
delay  line.  The  roots  of  these  weight  vectors  in  the  Doppler  domain  give  us 
an  estimate  of  the  Doppler  frequencies.  A  standard  polynomial  root  finding 
routine  is  utilized  for  this  purpose.  Knowing  the  frequencies,  the  MLM 
algorithm  is  then  applied  to  get  the  relative  power  levels  and  thus  estimate 
the  power  of  the  spectrum.  Once  again,  we  have  a  dramatic  reduction  in 
computer  processing  burden  through  having  an  estimate  of  the  particular 
Doppler  frequencies. 

The  final  step  consists  of  plotting  the  estimated  power  spectral  lines  at 
their  estimated  target  range  locations,  resulting  in  the  "digital"  3D  image  of 
Fig.  16(a).  This  method  overcomes  all  three  of  the  drawbacks  noted  for  the 
MLM/ LPF  Algorithm  images,  and  opens  up  the  possibility  for  practical 
technology  development. 

Triangle  Model  of  Ten  Point-Targets 

To  demonstrate  performance  against  a  more  complicated  rotating  object,  the 
final  simulations  consist  of  ten  point-targets  arranged  in  the  shape  of  a 
triangle  as  shown  in  Fig.  17.  TABLE  3  contains  the  various  parameter  values 
which  were  deliberately  chosen  to  result  in  close  target  spacings  and  to  fill 
three  adjacent  rangebins  with  combinations  of  three  or  four  targets  each. 
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Also,  note  that  five  of  the  targets  have  an  SNR  of  0  dB.  Data  samples  from 
264  prf  (total)  were  processed. 

Range  estimate  plots  are  given  in  Fig.  18  for  four  different  algorithms: 
Fourier,  MLM,  Wiener  SLC,  and  MUSIC  [30].  It  is  of  interest  to  note  the 
progression  of  resolution.  DFT  processing  gives  us  the  overall  "power 
envelope"  of  the  ten  targets  in  range;  MLM  processing  separates  them  into 
three  crude  groups;  Wiener  SLC  refines  the  estimate  further  by  indicating  that 
at  least  four  targets  are  present;  and  MUSIC  has  sufficient  sensitivity  to 
indicate  that  at  least  five  targets  are  present,  plus  a  number  of  extraneous 
artifacts.  The  latter  algorithm  is  a  sensitive  high-resolution  technique 
based  upon  the  noise  eigenvectors  of  the  covariance  matrix.  References 
[26,30]  are  recommended  for  those  readers  who  are  interested  in  further 
details  on  the  MUSIC  algorithm. 

Fig.  19  then  illustrates  the  two-dimensional  radar  images  processed  via  the 
same  three  algorithms  as  used  for  Fig.  16.  The  conventional  Fourier  transform 
plotted  in  Fig.  19(c)  again  serves  as  a  baseline  reference  to  illustrate  the 
results  of  conventional  processing.  Here  we  can  barely  discern  a  crude 
triangular  shape  within  which  at  least  six  components  are  evident.  The  MLM 
algorithm  plotted  in  Fig.  19(b)  successfully  resolves  all  ten  targets, 
although  accompanied  by  some  smearing,  location  displacements,  and  the 
"corrugation"  effect  refered  to  earlier.  Recognition  of  the  overall 
triangular  shape  and  component  scatterers  is  accomplished.  The  digital  image 
plotted  in  Fig.  19(a)  then  refines  the  estimate  and  enhances  the  component 
scatterers.  Note  that  the  digital  image  does  contain  a  few  extraneous 
artifacts,  but  these  are  not  considered  to  be  serious  defects  because  such 
artifacts  tend  to  move  around  (or  scintilate)  as  successive  images  are 
processed,  whereas  the  true  scatterers  remain  stable. 
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Triangle  Model  with  Coherence  Present 


In  section  4  we  alluded  to  the  deleterious  phenomenon  defined  as  "coherent 
scatterers”  wherein  two  or  more  closely-spaced  scatterers  have  the  same 
Doppler  frequency  as  seen  by  the  radar.  Our  final  simulation  example  is 
deliberately  designed  to  illustrate  such  a  case.  TABLE  4  contains  the  various 
parameter  values  for  the  ten  targets,  and  it  will  be  noted  that  coherence 
exists  (same  Doppler)  between  targets  No.  2  and  No.  10,  No.  4  and  No.  9,  and 
No.  6  and  No.  8.  Thus,  three  pairs  of  scatterers  within  the  triangle  are 
coherent . 

Fig.  20  illustrates  the  radar  image  results,  utilizing  the  same  three 
algorithms  as  before.  The  conventional  Fourier  transform  plotted  in  Fig. 
20(c)  gives  us  a  crude  triangular  shape  in  which  it  is  difficult  to  pick  out 
any  component  parts  beyond  the  corners  of  the  triangle.  The  MLM  algorithm 
plotted  in  Fig.  20(b)  does  a  little  better  in  resolving  component  parts,  but 
note  that  it  also  generates  some  distracting  outboard  artifacts  from  the 
processing.  The  digital  image  plotted  in  Fig.  20(a)  then  refines  the  estimate 
and  enhances  the  component  scatterers,  but  still  includes  the  rather  strong 
outboard  artifacts.  Careful  study  of  several  coherence  cases  has  revealed  the 
following  deleterious  effects: 

a)  the  positions  of  such  scatterers  may  be  shifted  in  both  range  and/or 
Doppler . 

b)  such  scatterers  may  blend  to  appear  as  a  single  scatterer  at 
intermediate  positions. 

c)  the  processed  power  levels  may  not  be  correct. 

d)  strong  artifacts  may  be  generated  at  range/ doppler  positions  where  no 
scatterer  exists. 
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Fig.  20  contains  ail  of  the  above  effects,  and  the  net  result  is  distortion 
of  the  processed  image.  To  overcome  these  effects,  it  is  necessary  to  employ 
a  decorrelation  method  such  as  the  carrier  shift  techniques  discussed  in 
Section  4. 

6.  Conclusions 

High-resolution  optimal  estimation  techniques  have  been  applied  to  the 
problem  area  of  radar  imaging  of  rotating  objects.  Typical  digital  range- 
Doppler  processor  operations  are  described,  utilizing  several  of  the  many 
spectral  analysis  estimation  techniques  available  in  the  literature.  Their 
application  to  the  range  domain  represents  new  R&D  work,  wherein  multiple 
scatterers  within  one  rangebin  have  been  resolved  without  increasing  the 
bandwidth.  Quality  ISAR  images  have  been  obtained  from  simulated  radar  data 
generated  from  point-target  models  of  rotating  objects.  These  simulation 
results  demonstrate  the  feasibility  for  achieving  superresolution  in  the 
Doppler  domain,  the  range  domain,  and  both  domains  jointly.  The  capability 
for  focusing  upon  a  small  section  of  the  overall  image  and  achieving  finer 
detail  (magnification)  is  particularly  useful.  It  is  concluded  that 
superresolution  techniques  offer  a  viable  complement  to  conventional  DFT-ISAR 
image  processing,  and  should  permit  either  higher  resolution  images  from  the 
same  data  samples,  or  equal-quality  images  from  significantly  fewer  data 
samples. 

Further  research  is  needed  in  several  areas  including  algorithm 
development,  algorithm  comparison,  meaningful  resolution  comparison  criteria, 
pulse-compression  waveform  design/comparison,  adaptive  "focussing"  upon 
dominant  scatterers,  techniques  for  reducing  computer  burden,  the  relationship 
between  image  quality  and  number  of  data  samples,  the  special  problem  area  of 


36 


coherent  target  scatterers,  rangebin  entering/ leaving  effects,  and  the 
evaluation  of  these  techniques  when  processing  experimental  radar  data 
samples . 
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TABLE  1 

ROTATING-BOOM  MODEL  for  SIMULATED  DATA 


fr=0.63 

rpm, 

prf=418  Rangebin=31. 

25  X  9= 

90  deg.  at 

Target 

No. 

SNR  Distance 

in  dB 

Doppler 
in  Hz. 

Doppler 
in  %prf 

Doppler 
in  deg. 

1 

0 

•  +45 . 0 

-185.6 

-44.4 

-159.8 

2 

0 

+30.0 

-123.7 

-29.6 

-106.5 

3 

0 

+19.0 

-78.3 

-18.7 

-67.5 

4 

23 

+10.0 

-41.2 

-9.9 

-35.5 

5 

23 

+9.2 

-37.9 

-9.1 

-32.7 

6 

23 

+1.6 

-6.5 

-1.6 

-5.6 

7 

23 

+0.8 

-3.3 

-0.8 

-2.8 

8 

23 

0.0 

0.0 

0.0 

0.0 

9 

0 

-5.1 

+20.9 

+5.0 

+18.0 

10 

20 

-6.8 

+27.8 

+6.7 

+24.0 

11 

10 

-8.4 

+34.8 

+8.3 

+30.0 

12 

20 

-10.1 

+41.8 

+10.0 

+36.0 

13 

0 

-19.0 

+78.3 

+18.7 

+67.5 

14 

0 

-30.0 

+123.7 

+29.6 

+106.5 

15 

0 

-45.0 

+185.6 

+44 . 4 

+159.8 
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TABLE  2 


ROTATING-BOOM  MODEL  for  SIMULATED  DATA 
fr=0.1  rpm,  prf=200,  Rangebin=150  A,  0=45  deg.  at  t=0 


4. 

Target 

SNR 

Distance 

Doppler 

Doppler 

Doppler 

No. 

in  dB 

d. 

i 

in  Hz. 

in  %prf 

in  deg. 

1 

10 

+1.0 

-3.78 

-1.89 

-6.81 

2 

10 

0.0 

-1.56 

-0.78 

-2.81 

3 

10 

-1.0 

+0.66 

+0.33 

+1.19 

Note : 

Doppler 

offset  of  -1. 

5625  Hz. , 

all  three 

targets 
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f  ■  0.1  rpm,  prf=  200,  Rangebin=  150  wavelengths,  t=  0  start  time, 
Doppler  offset  of  -1.5625  Hz.  all  ten  targets. 
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f  ■  0.1  rpm,  prf=  200,  Rangebin=  150  wavelengths,  t=  0  start  time, 
Doppler  offset  of  -1.5625  Hz.  all  ten  targets. 
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SHAPED  REFERENCE 
FILTER  OUTPUT 


e 

F 

i 

ii 

t 

E 

I 

E 

£ 

L 


FFE0UENC7  ID  Z-FLhNE  I'E'iREE> 

Fig.  4  —  Multiple  sub-bands  selected  from  DFT  which  connects  into  the  even-numbered  taps  along  the  64-tap  delay  line 
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ROTATING-BOOM 


FREQUENCY  IN  Z-»IANE  DEGREES 
Fig.  7  —  (a)  Conventional  Fourier  transform 


Fig.  7 


(c)  Sub-band  MLM  Algorithm. 

—  Typical  spectrum  estimates  for  nine  sinusoids  within  passband;  the  five  on  the  left  are  23  dB;  the  four  on  the  right 
are  0,  20,  10  and  20  dB  signaJ-to-noise  ratio.  140  data  samples  processed 
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12 


13 


14 


Fig.  8  —  “Merry-Go-Round”  of  point-targets  (24)  equally  spaced,  utilized  for  generating  simulated  radar  data  samples 
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3D  ISAR  images  of  “Merry-Go-Round"  routing  object,  processed  per  MLM  Algorithm,  32  up  delay  line. 

100  prf  processed 


:R  In  DB 


POWER  In  DB  POWER  in  DB 


(b)  MLM  Algorithm 

Fig.  13  —  Target  range  estimate  plots  obtained  from  two  high-resolution  algorithms,  264  prf  processed 
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RANGEBINS 


RANGE BINS 


<b)  MLM  Algorithm  outputs  averaged  over  264  prf. 

MLM  Algorithm  “Copy”  weights  filter  response,  as  applied  to  the  rangebin  region  from  12.0  to  18.0 
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Fig.  17  —  Rotating  object  consisting  of  ten  point-targets  arranged  in  the  shape  of  a  triangle 
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